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Abstract. Interfaces between two fluids are ubiquitous and of special 
importance for industrial applications, e.g., stabilisation of emulsions. 
The dynamics of fluid-fluid interfaces is difficult to study because these 
interfaces are usually deformable and their shapes are not known a 
priori. Since experiments do not provide access to all observables of 
interest, computer simulations pose attractive alternatives to gain in- 
sight into the physics of interfaces. In the present article, we restrict 
ourselves to systems with dimensions comparable to the lateral inter- 
face extensions. We provide a critical discussion of three numerical 
schemes coupled to the lattice Boltzmann method as a solver for the 
hydrodynamics of the problem: (a) the immersed boundary method for 
the simulation of vesicles and capsules, the Shan-Chen pseudopoten- 
tial approach for multi-component fluids in combination with (b) an 
additional advection-diffusion component for surfactant modelling and 
(c) a molecular dynamics algorithm for the simulation of nanoparticles 
acting as emulsifiers. 



1 Introduction 

Interfaces are ubiquitous in soft matter systems and appear in various shapes and 
sizes. Prominent examples are vesicles (closed membrane defined by a bilayer of phos- 
pholipid molecules [Tj), capsules (closed polymeric membrane [2J, and all kinds of 
biological cells [2j. In these cases, the interface is an additional material whose consti- 
tutive behavior has to be specified. For example, vesicle membranes are incompressible 
and viscous pQ whereas capsule membranes resist shear and can be compressible [3] . 
Understanding the dynamics of membranes is important for disease detection by mea- 
suring mechanical properties of living cell membranes |3] , targeted drug delivery [5] , 
and predicting the viscosity of biofluids, such as blood [Sj. Typical applications are 
lab-on-chip devices for particle identification and separation [7J. 
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Other classes of fluid-fluid interfaces can be found in emulsions (liquid drops sus- 
pended in another liquid), foams (gas bubbles separated by thin liquid films), and 
liquid aerosols (liquid drops in gas) where at least two immiscible fluid phases are 
mixed. The interface is then defined by the common boundaries of the phases. Emul- 
sions are of central importance for food processing (e.g., milk, salad dressings) [5], 
cosmetics and pharmaceutics (e.g., lotions, vaccines, disinfection) [5], and enhanced 
oil recovery [TU] , 

Emulsions and foams are usually unstable. Drops and bubbles tend to coalesce 
gradually, which reduces the interfacial free energy. Conversely, droplet breakup can 
be triggered by shearing the system and pumping energy into it. This is of practical 
importance for mechanical emulsification A significant amount of research has 
been focused on the question how to stabilize emulsions, thus avoiding or at least 
retarding coalescence. The traditional approach is to add surfactants to an emulsion 
[12j . Surfactants are amphiphilic molecules for which it is energetically favorable to 
accumulate at the interface, which in turn leads to a decrease of surface tension 
and prevents demixing. Under some circumstances, it is possible to find so-called 
mesophases where the phases percolate, penetrate each other, and are separated by 
a interfacial surfactant monolayer. These mesophases can be periodic such as the 
diamond or gyroid minimal surfaces |12H14j . 

An alternative route to emulsion stabilization is by using solid nanoparticles: these 
particles also accumulate at the interface where they replace segments of fluid-fluid 
interface by fluid-particle interfaces and thus lower the interfacial free energy. How- 
ever, the physical mechanism is different from that of stabilization by surfactants, 
and the surface tension is unchanged by the nanoparticles |15H17j . For nanoparticle- 
stabilized emulsions, one distinguishes between the so-called "Pickering emulsions" 
[T5J, [TS] and "bijels" (bicontinuous interfacially jammed emulsion gels) [501 I2J]- The 
former is an emulsion of discrete droplets in a continuous liquid. In the latter, both 
phases are continuously distributed. Besides their primary task to act as emulsifiers, 
the nanoparticles may additionally be designed as, for example, ferromagnetic |22| or 
Janus particles [23] . 

The macroscopic dynamics of complex fluids strongly depends on the microscopic 
properties of interfaces present in these systems. The reason is that the ratio between 
interface surface and bulk volume is usually large. The interfaces act as additional 
boundary conditions restricting the dynamics of the fluids. Quantities of interest are 
the interface's shear and dilatational viscosities and the surface dilatational modulus 
if the interface is compressible (which is — in contrast to bulk rheology — often the 
case). Interface models usually have to obey mass, momentum, and energy conserva- 
tion. If isothermal or athermal situations are sufficient, the energy conservation can 
be relaxed. For a review of experimental methods characterizing interfaces in soft 
matter and finding thcrmodynamically consistent constitutive laws, we refer to [24 . 

Fluid-fluid interface problems are hard to solve analytically since the interface is 
usually deformable and its shape not known a priori. Therefore, the interface dynam- 
ics is fully coupled to that of the ambient phases. Interface deformations are rarely 
small and linear, which results in complex and time-dependent boundary conditions. 
Analytical solutions are only available for academic cases (e.g., [25l[26]). These cir- 
cumstances call for numerical methods and computer simulations which can also 
provide access to observables not traceable in experiments such as local interface cur- 
vature or fluid and interface stresses. Therefore, computer simulations can be used to 
complement experiments. 

Engineering applications call for a better understanding of the link between the 
microscale and the emergent macroscopic behavior of these systems. For example, 
it is an open question how to generally predict the macroscopic behavior (such as 
the effective viscosity) of a complex fluid based on its microscopic properties (e.g., 
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(a) continuum (b) mesoscopic (c) discrete 



Fig. 1: Sketch of coarse-graining levels for an interface between two fluids (fluids 1 
and 2). In panel (a), the interface (green line) is explicitly tracked in a Lagrangian 
manner, and the interface properties have to be provided in form of a constitutive 
law. There is no direct interaction between the fluid species. The substructure of 
the interface is not known (continuum picture). For a surfactant-stabilised droplet as 
shown in panel (b), the interface is diffuse with an additional dynamical surfactant 
field (green). The local mutual interactions between the surfactant and the two fluid 
components have to be specified, but the microscopic surfactant properties are not 
required (mesoscopic picture). In panel (c), the nanoparticles (green) are explicitly 
resolved (discrete picture). Their individual motions are tracked, and their geometrical 
and interaction properties have to be specified (e.g., shape, contact angles). The 
interaction of the fluid species has to be defined as well. For (b) and (c), the choice 
of the mutual interactions leads to an emergent macroscopic interface behaviour. 



surface tension or dilation modulus of the interface, wettability of the nanoparticles) 
and how to manipulate these properties in order to tailor new materials with spe- 
cific macroscopic behavior [2"?H2l?] . What are the underlying physical mechanisms for 
emulsion stabilization due to surfactants and nanoparticles, and how can their costs 
and environmental sustainability be optimized? How can emulsions with a desired 
droplet size distribution be produced? How do drops, vesicles, and capsules behave in 
non-trivial flow geometries, and how can they be separated based on their interfacial 
properties? 

In this article, we focus on vesicles, capsules, and emulsions in systems with dimen- 
sions comparable to the lateral extension of the interfaces, e.g., a volume containing 
0(1 — 100) objects. At this scale, fluid dynamics is dominated by effects due to vis- 
cosity, surface tension or elasticity, and inertia. Gravity is usually not relevant in 
microfluidics because the capillary length for water droplets in air is typically of the 
order of a few millimeters. Due to the small system size, the total amount of energy 
dissipated is small and heat is almost instantaneously transported out of the system 
by heat conduction. It is therefore an excellent approximation to assume isother- 
mal systems with constant viscosity and surface tension. The relative strength of the 
three major contributions (viscosity, surface tension, inertia) can be described by di- 
mensionless parameters, such as the capillary number (Ca, viscous stress vs. surface 
tension or elasticity), the Reynolds number (Re, inertial vs. viscous stresses), or the 
Weber number (We, inertial stress vs. surface tension). Note that the three numbers 
obey We oc Re x Ca and are not independent. 

In the current article, the fluid phases are treated as Newtonian fluids simulated 
by the lattice Boltzmann method (detailed in section [2]). Depending on the physical 
conditions and the desired level of coarse-graining, it may be sufficient to treat the 
interface as an effective two-dimensional material obeying a well-defined constitutive 
law or it may be necessary to resolve the substructure of the interface (cf. fig. [I]). 
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We consider three cases: In the continuum description for vesicle- and capsule- mem- 
branes, the immersed boundary method is coupled to the single-component lattice 
Boltzmann method (section 2.1 1. The lattice Boltzmann Shan-Chen multi-component 
model is the basis for the remaining two cases. In th e first, a surfactant concentration 
field is employed (mesoscopic approach, section 2.2 1, in the second, explicitly resolved 
solid nanoparticles interact with the fluid phases (discrete approach, section 2.3). In 
section |3j some exemplary results from simulations of a vesicle under shear now, 
separating red blood cells and platelets, surfactant stabilized interfaces, and particle 
stabilized interfaces are presented. Finally, we conclude in section [4] 



2 Numerical approaches 

The physical problem consists in capturing the dynamics of an interface separating 
distinct fluids. Most of the macroscopic properties of soft matter systems, for exam- 
ple the rheology of an emulsion, result from the reorganisation and structuring of 
interfaces at the microscale. In the simplest situation, there are two different fluids 
separated by an interface. Mathematically, this corresponds to two domains separated 
by a free-moving boundary. There are two classes of numerical approaches to track the 
motion of an interface: (i) either front tracking (e.g., the boundary integral method 
or the immersed boundary method) in which the motion of marker points attached 
to the interface is tracked, or (ii) front capturing (e.g., the phase field method or 
level set methods) where a scalar field is used as an order parameter, indicating the 
composition of the fluid at a given point. For example, the order parameter may take 
the values and 1 in the bulk regions of fluid component A and B, respectively, and 
it is smoothly varying in the intermediate regions. The interface is located where the 
scalar field assumes a well-defined value, e.g., 0.5 in the above example. The time evo- 
lution of the order parameter is governed by an advection-diffusion equation coupled 
to the fluid flow. In both approaches (front tracking and capturing), the momentum 
of the incompressible Newtonian fluid obeys the Navier-Stokes (NS) equations, 

p + u ■ Vuj = -Vp + r?V 2 u, V u = 0. (1) 

The mass density and dynamic viscosity of the fluid are denoted by p and r\\ u and 
p are its velocity and pressure fields and t denotes the time. For the solution of the 
NS equations, boundary conditions at the interface have to be considered. However, 
the location of the interface is on its own an additional unknown quantity. Its motion 
results from its hydrodynamical interaction with the surrounding fluid. Therefore, to 
obtain the dynamics of the interface, information about the flow is required. In the 
present article, instead of solving the NS equations directly, we alternatively use a 
mesoscopic approach: the lattice Boltzmann method (LBM). 

The LBM has gained popularity among scientists and engineers because of its 
relatively straightforward implementation compared to other approaches such as the 
finite element method (for reviews, see [3"0Tf3"2"] L In the LBM, the fluid is considered as 
a cluster of pseudo-particles that move on a lattice under the action of external forces. 
To each pseudo-particle is associated a distribution function the main quantity 
in the LBM. It gives the probability to find a pseudo-particle at a position r with 
a velocity in direction e;. In the LBM, both the position and velocity spaces are 
discrctiscd: Ax is the grid spacing, and e$ are the discretised velocity directions. 
There arc different types of lattices available: For 2D simulations, we use the so-called 
D2Q9 lattice with nine velocity directions (see also Fig. [2]ja)); for 3D simulations, the 
D3Q19 lattice with nineteen velocity directions is employed [35]. The time evolution 
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of fi is governed by the so-called lattice Boltzmann equation, 

fi(r + ei At, t + At)- Mr, t) = fj t , (2) 

where At is the discrete time step. The lattice Boltzmann equation ^ consists of two 
parts: (i) the advection part (left-hand side) and (ii) the collision part (right-hand 
side) with collision operator f2i specifying the collision rate between the fluid pseudo- 
particles. It can be approximated by the Bhatnagar-Gross-Krook (BGK) operator, 
fli = — (fi — / Z eq )/T, which describes the relaxation of fi towards its local equilibrium, 
/, eq , on a time scale r. The relaxation time is related to the macroscopic dynamical 

viscosity 77 via rj = pc^ (r — |) , where c s = 1 /V% is the lattice speed of sound. The 
equilibrium distribution is given by a truncated Maxwell-Boltzmann distribution, 



fP( P ,u)=LO lP 



c, • u (c, • u) 2 u • u 

1 + — h — — 

c 2 2 C 4 a 



(3) 



where the Wi are weight factors resulting from the velocity space discretisation. The 
hydrodynamical macroscopic quantities are computed using the first and the second 
moments of f^. (i) the local pressure p = pc% = c 2 , fi ancl (ii) the local velocity 
u = fi&i/ p. Even though a conversion of the units to SI units is straightforward, 
it is common practise to report all quantities in lattice units. Non-dimensionalisation 
is possible by choosing suitable combinations of Ax, At, and the fluid density. 

The LBM can also handle multi-phase and multi-component fluids and a number 
of corresponding extensions of the method have been published in the past [34H39] , 
In the Shan-Chen multi-component model [33], a system containing N miscible or 
immiscible fluids (with index a) are described by N sets of distribution functions 
/cr.j, one for each species. As a consequence, N lattice Boltzmann equations with 
relaxation times r a have to be considered. The interaction between the fluid species 
is mediated via a local force density, 

F ff (r,t) = -& a (r,t)J2 W$> CT '(r',«)(r'-r). (4) 

a' r' 

Here, the contributions of all species a' arc taken into account, and the force density 
acting on species a is obtained. The sum runs over all neighbouring sites r' of site 
r. The factor g aa i is the coupling constant defining the interaction strength between 
species a and a'. It is related to the surface tension between both species. Self- 
interactions (a' = a) are also allowed. The pseudopotential ^ a is a function of the 
density p a , in our case ^ a = 1 — exp(— Pa ). Its shape is related to the equation of 
state. There are different ways to incorporate forces such as those in Eq. (4]) into the 
LBM. One possibility is the so-called velocity shift in the equilibrium distribution 
function: For the computation of the equilibrium distribution /° q - of component a, u 
in Eq. ^ is replaced by (u oq + T a F a /p a ), where p a is the density of component a 
and u eq = J2 a P " T u ' y I ^ 1S the common equilibrium velocity of all components. 
The flow velocity of the fluid is then given by u = 53<r(Si /o-,* e « + P<r/2)//0. The 
Shan-Chen model belongs to the class of front capturing methods. It is suitable to 
track interfaces for which the topology evolves in time, for example, the breakup of 
a droplet. 

In section [3] we give examples where the LBM is used as a fluid flow solver to 
study the dynamics of complex interfaces. 



2.1 Continuum interface models 



In this section, we present the immersed boundary method (IBM) as a front tracking 
approach [40j . Such methods are mostly employed when the interface is formed by an 
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Fig. 2: Two-dimensional illustrations of the lattice Boltzmann propagation, immersed 
boundary interpolation, and bounce-back, (a) During propagation, the populations 
fi move to their next neighbors (gray lattice sites and arrows) . Within the immersed 
boundary method, an interface node (here: node i at position rj(t)) is coupled to the 
single-component lattice fluid. The range of the interpolation stencil (here: 4 x 4) is 
denoted by the square region. All lattice sites within this region have to be considered 
during velocity interpolation and force spreading, (b) An ellipsoidal nanoparticle is 
located at a fluid-fluid interface. One fluid component is indicated by red, the other 
by blue circles. The colour gradient illustrates the interface region. Squares denote 
lattice sites inside the particle. Populations are bounced back at points (shown as 
small circles) located half-way between neighbouring fluid and particle sites. This 
gives rise to an effective staircase description of the particle shape. 



additional continuous material whose constitutive behaviour (e.g., elasticity, viscos- 
ity) is assumed to be known. Examples are capsules, vesicles, or biological cells. It is 
usually not necessary to resolve the physical mechanisms which lead to this consti- 
tutive behaviour since there is a significant scale separation (if one is only interested 
in the mechanical properties on the macroscale) . The IBM and related front-tracking 
approaches can also be used to describe the motion of the interface even when it is 
not formed by an additional material. For example, the interface between a gas and 
a liquid phase may be tracked [IT] , 

Within the IBM, the interface is considered sharp (zero thickness) and is repre- 
sented by a cluster of marker points (nodes) which constitute a moving Lagrangian 
mesh (cf. Fig. [2ja)). This mesh is immersed in a fixed Eulerian lattice representing 
the fluid. The question is how to predict the motion of the interface in time and how it 
affects the fluid. To consider correct dynamics, a bi-directional coupling of the lattice 
fluid and the moving Lagrangian mesh has to be taken into account. On the one hand, 
the interface is modelled as an impermeable structure obeying the no-slip condition at 
its surface. It is assumed that the flow field is continuous across the interface and that 
the interface is massless. Therefore, the interface is moving along with the ambient 
fluid velocity. On the other hand, a deformation of the interface generally leads to 
stresses reacting back onto the fluid via local forces. As an example, interfacial stresses 
are caused by local bending of a fluid-fluid interface (surface tension) or shearing of a 
capsule membrane (shear deformation). The stresses depend on the chosen constitu- 
tive behaviour of the interface and are not predicted within the IBM itself. This makes 
the IBM a flexible coupling scheme as it does not dictate specific material properties. 
The two-way coupling is accomplished in two main steps: (i) velocity interpolation 
and Lagrangian node advection and (ii) force spreading (reaction). 
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Interpolation: First, the fluid flow field is computed such that the velocity u is 
known at each lattice site x. In principle, any Navier-Stokes solver is applicable; in 
our case, we opted for the LBM, resulting in an immersed boundary-lattice Boltz- 
mann method. Afterwards, the fluid velocity is interpolated to obtain its value at the 
position Tj of each membrane node i which equals the velocity of that node: 



The interpolation stencil A(x — r^) is a suitable discretised approximation of Dirac's 
delta function [40]. For numerical efficiency, its range should be finite. Peskin sug- 
gested several stencils with different interpolation ranges, e.g., taking into account 
two, three, or four lattice sites along each coordinate axis [20]. While the 2-point 
stencil is numerically faster and leads to a sharper description of the interface, the 3- 
and 4-point stencils result in smoother interpolated velocity fields [22]. Finally, the 
interface is advected by updating the position of each membrane node using an Euler 
scheme, r<(i + At) = r^t) + ii(t)At. 

Reaction: By advecting the interface shape, it is generally deformed. The new 
shape of the interface is not necessarily its equilibrium shape which would minimise 
its energy. Therefore, each node exerts a reaction force Fj (which is computed from the 
known deformation state and the constitutive interface properties) on its surrounding 
fluid. This force is distributed (also called "spread") to the fluid using the same stencil 
as for the velocity interpolation, 



This force has to be taken into account by the Navier-Stokes solver as a local accel- 
eration in the next time step. 

2.2 Mesoscopic interface models 

Surfactants have been treated numerically for many years and as such, the models and 
descriptions have evolved over time. Even though amphiphilic surfactant molecules 
derive many important properties from their structure at microscopic scales (the pres- 
ence of molecular hydrophilic "head" and hydrophobic "tail" groups) , modelling them 
as individual molecules is generally unfeasible, computationally, if one is interested 
in effects on macroscopic or mesoscopic scales. For example, in the case of surface 
tension reduction at a fluid-fluid interface, this would not supply any detail that is 
not obscured by the rheology; thus, surfactants are often modelled with some de- 
gree of coarse-graining. This makes this strategy well-suited to be coupled to other 
mesoscopic methods. 

Simulations based on lattice-gas-automata (LGA) have been used since the nine- 
ties. As a strongly simplified scheme, only the effect that the surfactant has on surface 
tension can be modelled |43j . Gompper and Schick have introduced a model which 
allows the surfactant to have orientational degrees of freedom, as well as translational 
ones [14| , which due to the shape and nature of the amphiphiles is a necessary prop- 
erty to simulate. Just as the LBM has evolved from LGA, the attendant surfactant 
extensions have found their way into LB as well. These implementations, built as an 
extension of LB, are — unlike LGA-based models — not suitable for simulating phe- 
nomena that strongly depend on effects caused by discrete particles. However, they 
are perfectly adequate in the area of bulk hydrodynamic phenomena which is of inter- 
est here. For example, the model developed by Chen et al. [24H2B] uses self-consistent 
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forces to describe the interactions. This method introduces the surfactant as both 
an additional scalar field to model the local densities and a vector field to model 
local average dipole orientations (which can vary continuously). These fields are two- 
way coupled to the non-amphiphilic species in the system through Shan-Chen-type 
interaction forces. A similar force also describes a self-interaction to model the at- 
traction between two amphiphile tails and repulsion between a head and a tail. This 
method is rather simple in its implementation, and as it assumes only interactions 
with neighbouring lattice sites, it is local and fast. However, the effective surface 
tension reduction that can be effected through this model is rather modest [17]. 

Alternative methods built upon the LBM include the Ginzburg-Landau free energy- 
based model with two scalar order parameters introduced by Lamura et al. [471 and 
an extension to Chen's method which includes mid-range interactions by using more 
than one Brillouin zone in the calculations [IS] 25] ■ This increases the range of sur- 
face tension reductions that can be achieved at the cost of reduced computational 
efficiency. 



2.3 Discrete interface models 

If there is no clear scale separation between immersed particles and the lateral inter- 
face extension, it may be necessary to model the particles explicitly. Here, we consider 
an ensemble of particles with well-defined wetting behaviour. The particles themselves 
may have a complex shape and can be provided with a constitutive model (of which 
the simplest is a rigid particle model as we use it here) . There are several approaches 
to simulate a system of two immiscible fluids and particles. One example which has 
recently been applied by several groups is Molecular Dynamics (MD) coupled to the 
LBM 17, 20, 22, 50-53.. The advantage of this combination is the possibility to resolve 
the particles as well as both fluids in such a way that all relevant hydrodynamical 
properties are included. The particles are generally assumed to be rigid and can have 
arbitrary shapes, where we restrict ourselves to spheres and ellipsoids. 

For the fluid-particle coupling, the particles are discretised on the lattice: sites 
which are occupied by a particle are marked as solid. Following the approach pro- 
posed by Ladd [54] , populations propagating from fluid to particle sites are bounced 
back in the direction they came from (cf. Fig.^b)). In this process, the populations 
receive additional momentum due to the motion of the particles. In order to satisfy 
the local conservation laws, linear and angular momentum contributions are assigned 
to the corresponding particles as well. These in turn are used to update the particle 
positions and orientations. In the examples presented in this paper, a leapfrog inte- 
grator is used to solve the equations for linear and angular particle momenta. In the 
context of the Shan-Chen multi-component algorithm coupled to the particle solver, 
the outermost sites covered by a particle are filled with a virtual fluid corresponding 
to a suitable average of the surrounding unoccupied sites. This approach provides 
accurate dynamics of the two-component fluid near the particle surface. The wetting 
properties of the particle surface can be controlled by shifting the local density dif- 
ference of both fluid species by a given amount Ap. It can be shown that the contact 
angle 6 of the fluid-fluid interface at the particle surface has a linear dependence on 
Ap [50] [51]. Occasionally, when the particles move, lattice sites change from particle 
to fluid state, which has to be treated properly [T71 [5U] . 

If particles approach each other so closely that the lattice resolution is not suffi- 
cient to resolve their hydrodynamic interaction, an additional lubrication correction 
is required [5011511 [55] . Further, a Hertzian potential is used to approximate an elastic 
interaction between pairs of particles. For two overlapping spheres of radius R with 
mutual centre distance < 2R, it is defined as </>h = K U (2R - r l3 f' 2 [55] where 
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the force magnitude is controlled by K^l- A generalisation for ellipsoids, based on the 
method of Berne and Pechukas [5 7) , is explained in [51] . 

3 Case studies 

3.1 Continuum interface examples 

In this section, we show two exemplary applications of the immersed boundary-lattice 
Boltzmann method as introduced in section |2.1| First, we present two-dimensional 
simulations of the dynamics of a viscous vesicle in simple shear flow. In section [3. 1.2[ 
three-dimensional simulations of the separation of blood components are shown. 

3.1.1 Tumbling and tank-treading of viscous vesicles 

A vesicle is a closed phospholipid membrane filled with a fluid. Vesicles can be used 
to micro-encapsulate active materials for drug delivery. They also pose a biomimetic 
model for biological cells. One of the questions that attracted the interest of scientists 
in the past decades was: how does a vesicle behave dynamically under flow? There 
exists a rich literature related to this problem (e.g., [58 -60 ). Today it is known that a 
vesicle undergoes either tank-treading (steady motion) or tumbling (unsteady motion) 
when it is subjected to simple shear flow. In the tank-treading mode, the vesicle's main 
axis assumes a steady inclination angle 9 with the flow direction while its membrane 
undergoes a tank-treading like motion. In the tumbling mode, the vesicle rotates as 
a solid, elongated particle. 

Here, we reproduce these dynamical states via the combined immersed boundary- 
lattice Boltzmann method in two dimensions at negligible Reynolds numbers [5_T1 [53] . 
The vesicle encloses an internal fluid and is suspended in an external fluid. Both fluids 
are considered incompressible and Newtonian. In the present case, we are interested 
in viscous vesicles where the internal and external fluid viscosities differ in general. 
Their ratio A := 7 ^ L is called viscosity contrast. Numerically, we distinguish between 
the two fluids by the even- odd rule, an algorithm deciding whether a fluid site is 
inside or outside the vesicle [55]. Depending on the location, the local fluid viscosity 
is set to either n ln t or 77 oxt via two different lattice Boltzmann relaxation times (ri n t 
and T ext ). The vesicle membrane (interface) influences the ambient fluids (inside and 
outside) by exerting a local force on them. This force results from the resistance of 
the membrane to bending and stretching/compression |63j . It is computed on the 
membrane according to 



where n and t are the normal and tangential unit vectors, respectively. The membrane 
is characterised by the membrane bending modulus Kb and the local effective tension 
Q. The local curvature is c, and s is the arclength coordinate. The 4-point immersed 
boundary stencil is used to couple the dynamics of the membrane on the Lagrangian 
mesh and the fluid on the Eulerian lattice. 

Fig. [3] shows the dynamics of a vesicle in a simple shear flow with shear rate 
7. For a smaller viscosity contrast {A = 1), the vesicle tank-treads (Fig. [3]ja)). The 
tank-treading motion of the vesicle membrane generates a rotational flow inside the 
vesicle. By increasing the viscosity contrast above a certain threshold, we induce a 
transition from tank-treading to tumbling. Fig. p5lb) shows the tumbling motion of a 
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Fig. 3: Snapshots showing the dynamics of a viscous vesicle under shear flow, (a) 
Tank-treading: the vesicle assumes a steady inclination angle with the flow direction 
while its membrane undergoes a tank-treading motion, (b) Tumbling: the vesicle 
rotates as a rigid elongated particle. The six snapshots (from left to right) illustrate 
the time evolution of the inclination angle, cf. Fig. |4j 
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Fig. 4: Time evolution of the vesicle inclination angle 9 during the tank-treading 
(viscosity contrast A = 1) and tumbling (A = 15) modes. The time is scaled by the 
shear rate j. Some snapshots of the vesicle are shown in Fig. [3| 



vesicle with A = 15. In Fig. [4j we report the time evolution of the inclination angle 
9 for both cases. Initially, the vesicle main axis is parallel to the flow axis (9 = 
at t = 0). In the tank-treading mode, 9 increases until it assumes a steady value. 
Contrarily, in the tumbling mode, 9 varies periodically. 

It can be seen that the relatively simple immersed boundary-lattice Boltzmann 
method can be successfully employed to study vesicle dynamics in external flow fields. 
In particular, no remeshing of any kind is required. 



3.1.2 Flow-induced separation of red blood cells and platelets 

Efficient separation of biological cells plays a major role in present-day medical ap- 
plications, for example, the detection of diseased cells or the enrichment of rare cells. 
While active cell sorting relies on external forces acting on tagged particles, the idea 
of passive sorting is to take advantage of hydrodynamic effects in combination with 
the cells' membrane (interface) properties. It is known that deformable and rigid par- 
ticles behave differently when exposed to external flow fields. For example, in Stokes 
flow (Re = 0), deformable capsules (Ca > 0) show a strong tendency to migrate 
towards the centreline of a pressure-gradient driven channel flow. Therefore, it has 
been proposed to use the particle deformability as intrinsic marker to separate rigid 
and deformable particles in specifically designed microfluidic devices. We show that 
the separation of red blood cells (RBCs) and platelets, which has been shown recently 
via experiments |64j , can be simulated by means of the LBM and immersed boundary 
method. 
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(a) initial state (b) final state, Ca = 0.015 (c) final state, Ca = 0.3 

Fig. 5: Simulation snapshots of the separation of red blood cells (RBCs) and platelets. 

(a) Initially, the random suspension of 40 RBCs and 30 platelets is located in the 
bottom 30% of the channel (50 |xm diameter, denoted by two black lines). RBCs 
(red) are shown with reduced opacity to reveal the platelet positions (yellow) . Panels 

(b) and (c) show the lateral distribution of the cells after they have been moving 
downstream (rightwards) by about 20 mm on average. For rigid and tumbling RBCs, 
lateral separation of RBCs and platelets is not pronounced (b). The separation is 
more efficient when the RBCs are strongly deformable and are tank-treading in the 
vicinity of the walls (c). The corresponding volume fraction profiles are shown in Fig. 

E 

The RBCs and platelets are modelled as closed two-dimensional membranes with 
identical fluids in the interior and exterior regions [A = 1). In the present approach, 
a finite element method is used to compute the local surface shear deformation and 
area dilatation and the resulting forces [42 , 65] . The 2-point stencil has been used for 
IBM interpolation and spreading. 40 RBCs and 30 platelets are randomly distributed 
in the bottom 30% of a channel segment with 50 x 50 x 50 u\m 3 volume, cf . Fig. [5| The 
platelets can be considered as rigid ellipsoidal particles. The large RBC and platelet 
diameters are 16 and 5.5 lattice constants, respectively. The channel is bounded by 
two parallel walls and is periodic in the other two directions. A force / in rightward 
direction is used to drive the flow, resulting in a Poiseuille-like velocity profile. The 
average shear rate is defined as 7 := 2u/H where u = fH 2 /(8rj) is the central velocity 
in the absence of particles, H is the channel diameter, and ij is the external fluid 
viscosity. The RBC membrane (interface) is characterised by the capillary number 
Ca = rjjr/Ks. Here, r = 4u\m is the large RBC radius and n s is its shear elasticity 
which is about 5 |xNm _1 for healthy RBCs. Two simulations have been run, one with 
the normal elasticity, leading to Ca = 0.3 for the selected viscosity and driving force. 
In the other case, the RBC shear modulus has been increased by a factor of 20, and 
the capillary number was 0.015. This situation is typical for diseased RBCs (e.g., due 
to malaria or sickle cell anaemia). 

The initial simulation states are shown in Fig. [5^a) and Fig.[6|a), while the corre- 
sponding panels (b) and (c) show the state after about 20 mm downstream motion for 
different capillary numbers (0.015 and 0.3). It can be seen that RBCs and platelets can 
be efficiently separated when the RBCs are sufficiently deformable (Fig.[5](c) and Fig. 
(6jc)). The platelets marginate into the gap forming between the bottom wall and the 
RBC bulk. Contrarily, for the rigid RBCs, separation is visible but not pronounced 
(Fig. (5{b) and Fig. [f|b)). Additionally, the bulk of the RBCs moves faster towards 
the centreplane when the particles are more deformable, i.e., when Ca is larger. For 
Ca = 0.30, RBCs near the wall are observed to tank-tread; for Ca = 0.015, all are 
tumbling. The lateral migration timescale is several orders of magnitude larger than 
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Fig. 6: Local volume fraction <f> (averaged over flow and vorticity directions) of red 
blood cells (RBCs, dashed red lines) and platelets (solid black lines) as function of 
lateral position. The two walls are located at and H . The platelet volume fractions 
are multiplied by a factor of 10 for better visibility. The plots are rotated by 90° to 
enable a direct comparison of the profiles, (a) Initial state; state after about 20 mm 
downstream motion for (b) rigid and (c) soft RBCs. Although platelet margination 
away from the RBC bulk can be seen in both cases, the platelet separation towards 
the bottom wall is more efficient for softer RBCs. The three panels correspond to the 
snapshots in Fig. [5] 



the timescale for advection along the channel: the bulk of the soft RBCs has moved 
by only 15 u.m towards the centreplane after travelling about 20 mm in downstream 
direction. 

This case study is a striking example for the effect of interface properties on the 
overall flow behaviour of suspensions. It also shows that the immersed boundary- 
lattice Boltzmann method can be applied to suspensions of soft particles at moderate 
volume fractions. 



3.2 Mesoscopic interface examples 

A mesoscopic description of interfaces stabilised by amphiphilic molecules involves a 
systematic coarse-graining of the molecular degrees of freedom, while keeping the ef- 
fective mediated interactions between other fluids in the system. This can be obtained 
using the model described in section [2.2[ where the surfactant is simulated as a fluid 
phase together with a local mean-field dipolar force mediating the amphiphilic inter- 
actions on a mesoscopic level. The applicability of the method is being demonstrated 
in this section by two examples: first, we show how phase separation and emulsion 
stabilisation can be simulated using a ternary LB model. Second, the model is applied 
to study a self-organised amphiphilic mesophase, namely the gyroid mesophase. 



3.2.1 Phase separation and emulsification of amphiphilic mixtures 

In section [2. 2[ several mesoscopic methods to simulate amphiphilic fluids were men- 
tioned. Such fluids contain at least one species made of surfactant molecules and are 
in general complex fluids. By adjusting temperature, fluid composition or pressure, 
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Fig. 7: Volume rendered fluid densities for surfactant densities p s =0.00, 0.15, 0.30 
(from left to right). Only one of the immiscible fluid species is shown and different 
colours denote the interface and areas of high density of the visualised fluid. The 
surfactant (not shown) is mostly aligned at the interfaces and the second immiscible 
fluid component fills the void space. After 3 • 10 4 time steps the phases have separated 
to a large extent if no surfactant is present (left). Adding a small amount of surfactant 
(centre) causes the domain growth to slow down. For high concentrations (right) the 
growth process arrests and a stable bicontinuous microemulsion is formed (from [TTj). 



the amphiphiles can self-assemble and force the fluid mixture into a number of equi- 
librium structures or mesophases. These include lamellae and hexagonally packed 
cylinders, micellar, primitive, diamond, or gyroid cubic mesophases as well as sponge 
phases . The latter are generally not well ordered, but might have a well defined 
average size of fluid domains. In the context of our simulations they are called bicon- 
tinuous microemulsions since they are formed by the amphiphilic stabilisation of a 
phase-separating binary mixture, where the two immiscible fluid constituents occur 
in equal proportions. Here, the oil and water phases interpenetrate and percolate and 
are separated by a monolayer of surfactant at the interface. 

It has been shown by Langevin, molecular dynamics, lattice gas, and lattice Boltz- 
mann simulations that the temporal growth law for the size of oil and water domains 
in a system without amphiphiles follows a power law t 1 , with 7 being between 1/3 
and 1 |67j . Adding amphiphiles to a binary mixture of otherwise immiscible fluids 
which is undergoing phase separation can cause the separation to slow down. The 
domain growth behaviour then crosses over to a logarithmic growth (lni) e . A further 
increase of the surfactant concentration can lead to growth which is well described 
by a stretched exponential form A — B exp(— Ct s ). Capital and Greek letters denote 
fitting parameters and t is the time [S3^70|. 

In the simulations shown here we model two immiscible LB fluids and one surfac- 
tant species. We avoid effects due to variations of the fluid densities by keeping the 
total density in the system constant at 1.6 (in lattice units), while varying the sur- 
factant densities p s between 0.00 and 0.70. To depict the influence of the surfactant 
density on the phase separation process, Fig. [7] shows three volume rendered 256 3 
systems at surfactant densities 0.00 (left), 0.15 (centre), and 0.30 (right). After 3- 10 4 
time steps the phases are almost separated if no surfactant is present (left). Longer 
simulation times would result in two perfectly separated phases. If one adds some sur- 
factant (p s = 0.15, centre), the domains grow more slowly, visualised by the smaller 
structures in the volume rendered image. For sufficiently high amphiphile concentra- 
tions (p s = 0.30, right) the growth process arrests leading to a stable bicontinuous 
microemulsion with small individual domains formed by the two immiscible fluids. 
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Fig. 8: Linear (a) and logarithmic (b) representation of the average domain size L(t) 
for various surfactant densities p 8 (data reprinted from |71j). 



To investigate the influence of surfactant more quantitatively, we define the time 
dependent lateral domain size L(t) along direction i — x,y,z as 

where {kf(t)) = (J^ k fcf 5(k, t)) / Ek^C 4 '*)) ^ s tne secon d order moment of the 
three-dimensional structure function 5(k, t) = h 1 (i) | 2 with respect to the Carte- 
sian component i, () denotes the average in Fourier space, weighted by S(k.,t) and V 
is the number of sites of the lattice, </4(i) the Fourier transform of the fluctuations 
of the order parameter <j>' = <fi — (0), and ki is the ith component of the wave vector. 

In Fig. [8] the time dependent lateral domain size L(t) is shown for a number 
of surfactant densities p s between 0.00 and 0.50. Fig. ^a.) and (b) show identical 
data, but different scalings of the axes. In (a), we plot the data linearly in order to 
give a better impression of the time dependence of the growth dynamics. However, 
to determine the actual growth law, we provide a log-scale plot of the same data 
in Fig. ^h). For the first few hundred time steps, the randomly distributed fluid 
densities of the initial system configuration cause a spontaneous formation of small 
domains resulting in a steep increase of L(t). For p s — 0.00, domain growth does 
not come to an end until the domains span the full system. By adding surfactant we 
can slow down the growth process and for high surfactant densities, p s > 0.25, the 
domain growth stops after only a few thousand simulation time steps. By adding even 
more surfactant, the final average domain size decreases further. We fit our numerical 
data with the corresponding growth laws and find that for p s smaller than 0.15 L(t) 
is best fit by a function proportional to t a . For p s being 0.15 or 0.20, a logarithmic 
behaviour proportional to (\nt) e is observed. Increasing p s further results in L(t) 
being best described by a stretched exponential. These results correspond well with 
the findings in [70] . 



3.2.2 The cubic gyroid mesophase 

The ternary amphiphilic lattice Boltzmann method used here has been very suc- 
cessfully applied to study the dynamical self-assembly of a particular amphiphile 
mesophase, the gyroid [HI [T31 [551 EH [ZD ■ This mesophase is observed to form from 
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a homogeneous mixture of fluids, without any external constraints imposed to let the 
gyroid geometry emerge. This geometry is an effect of the mesoscopic fluid parameters 
only, which is remarkable, since it allows insight into the dynamics of the mesophase 
formation — in contrast to most other treatments to date which have focused on static 
properties or a mathematical description of the static equilibrium state |73l I74j . In 
addition to its biological importance, there have been recent attempts to use self- 
assembling gyroids to construct nanoporous materials [75] . The gyroid structure can 
be understood as two labyrinths mainly consisting of water and oil counterparts which 
are enclosed by the gyroid minimal surface at which the surfactant molecules accu- 
mulate. Fig. [9] shows a snapshot from a typical simulation of a gyroid mesophase 
and depicts the characteristic triple junctions defined by the volumes containing the 
majority of one of the individual fluid species. Multiple gyroid domains have formed 
and the close-up shows the extremely regular, crystalline, gyroid structure within a 
domain. 

In small-scale simulations, the gyroid mesophase will evolve to perfectly fill the 
simulated region, without defects. However, this can be accounted to finite size effects. 
As the lattice size is being increased, numerous small, well separated gyroid-phase re- 
gions or domains may start to form during the gyroid self-assembly. While these 
domains grow independently, they will in general not be identical, and might differ in 
orientation, position, or unit cell size. Thus, grain-boundary defects arise between gy- 
roid domains. Inside a domain, dislocations, or line defects might occur, corresponding 
to the termination of a plane of unit cells. There may also be localised non-gyroid 
regions, corresponding to defects due to contamination or inhomogeneities in the ini- 
tial conditions. Even with the longest possible simulation times, it is impossible to 
generate a 'perfect' crystal. Instead, either differently orientated domains can still be 
found or individual defects are still diffusing around. Understanding such defects is 
important for the comprehension of how these phases behave dynamically and how 
one could produce and utilise them experimentally [13] . 

Of further interest is the rheological behaviour of the gyroid mesophase. It has 
been shown that non-trivial rheological properties emerge when the gyroid phase is 
placed in an oscillatory shear flow with frequency uj /2tt. These include shear-thinning 
and even linear viscoelastic effects. Furthermore, the ternary amphiphilic LB model 
correctly predicts the theoretical limits for the moduli G'(ui), G"{uS) as u> goes to 
as well as a crossover in G'(oj), G"(to) at higher to [SHIEE]. It is remarkable that the 
model does not require any assumptions at the macroscopic level in order to predict 
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the formation and rheological behaviour of cubic mesophases. Instead, it provides a 
purely kinetic approach to the description of complex fluids. 



3.3 Discrete interface examples 

We now consider the effect of massive nanoparticles adsorped to a fluid-fluid interface 
as an example of the coupled lattice Boltzmann-molecular dynamics method explained 
in section |2.3| The first example shows the ordering of ellipsoidal particles at the 
interface of a spherical fluid droplet and compares it with the corresponding results 
of particles at a flat interface. This is followed by a case study highlighting how the 
presence of spherical nanoparticles affects the properties of a droplet in shear flow. 



3.3.1 Ordering of ellipsoidal particles at a fluid-fluid interface 

Particles adsorbed at a fluid-fluid interface can stabilise emulsions of immiscible fluids 
(e.g., Pickering emulsions where droplets of one fluid are immersed in another fluid 
and are stabilised by colloidal particles) [SSI HD • These colloids are not necessarily 
spherical such as a clay particle which has a flat shape. A simple approximation of 
such a particle shape is an ellipsoid. The adsorption of a single ellipsoid at a flat 
interface is the first step to understand the dynamics of such emulsions [ST] . The next 
step is the study of the behaviour of many ellipsoidal particles at an interface. In this 
section, the behaviour of an ensemble of elongated ellipsoidal particles with aspect 
ratio 77i = 2 at a single droplet interface is discussed [77] . The results are compared to 
the corresponding case of a flat interface. The initial particle orientation is such that 
the long axis is orthogonal to the local interface (cf. Fig. [lO"|a)). As discussed in [5T] . 
for the case of a single prolate ellipsoidal particle at a flat interface, the equilibrium 
configuration is a particle orientation parallel to the interface. We vary the coverage 
fraction \ = NA P /Ad which is defined as the ratio of the initial interface area covered 
by the N particles orientated orthogonal to the interface and the total surface area 
of the spherical droplet Ad- A p is the area which is covered by a single particle. 
Note that the covered area increases when the particles rotate towards the interface. 
The simulations discussed in this section have been performed for \ = 0.153 and 
X = 0.305. 



Fig. 10 'a) and (b) show the particle laden droplet almost at the beginning of 
the simulation (after 10 3 time steps) where all particles are perpendicular to the 
interface for both values of \- For the case of % = 0.305, the particles are comparably 
close to each other so that capillary interactions lead to the clustering motion of 
particles. The reason for these capillary interactions is a deformation of the interface 
caused by particles rotating towards the interface. In addition the rotating particles 
cause a small deviation of the droplet from its originally spherical shape which also 



leads to capillary interactions [75]. In Fig. 10{c), at the end of the simulation, the 
particles have completed a rotation towards the interface and stabilise the largest 
interfacial area possible. However, their local dynamics and temporal re-ordering is 
a continuous interplay between capillary interactions, hydrodynamic waves due to 
interface deformations and particle collisions. Inspired from liquid crystal analysis, 
we use the uniaxial order parameter S to characterise the orientational ordering of 
anisotropic particles, 

S = I (3 cos 2 0- lV (9) 



where 7? is the angle between the particle main axis and the normal of the interface 
at the point of the particle centre. The angular brackets denote an instantaneous 
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(a) t = 10 3 , x = 0.153 (b) t = 10 3 , x = 0.305 (c) t = 2 • 10 4 , X = 0.305 



Fig. 10: Snapshots of the particle-laden droplet for different times (in time steps) and 
coverage fractions x- One observes a higher global spatial order for the more dilute 
system in (a) as compared to the denser case in (b). 



ensemble average. An ideal spherical shape is assumed for the calculation of S. Fig. 



11 a) shows the time evolution of the order parameter S for particles at a droplet 
interface for \ — 0.153 and % — 0.305 and the corresponding cases for particles at a flat 
interface. Initially, the order parameter assumes the value S± = 1 which corresponds 
to perfect alignment perpendicularly to the interface. In the case of flat interfaces, 
the order parameter reaches the final value of S » Sn = —0.5. This corresponds to 
the state where all particles are aligned with the interface plane. For particle laden 



droplets (cf. Fig. 10 c)), a value of S s» —0.4 > S\\ is found, which indicates that 
the particles are not entirely aligned with the interface. Furthermore, S fluctuates in 
time. One possible reason for this deviation of S from S\\ is that the droplet does 
not have the exact spherical shape which is assumed for the calculation of S. This 
would lead to an overestimated value of S. The fluctuation of S can be explained by 
distortions of the droplet shape due to the dynamics of the particle rotations. Another 
difference between the curved and the flat interface is the time scale for the particle 
flip. However, one cannot see a significant effect of \ on the time development of the 
order parameter. 

To measure the spatial short range ordering of the particles, we define the pair 
correlation function 

^)=2^vEf?^-^i)^ (10) 

with discrete values for the distance r between particle centres. The number of par- 
ticles in the system is given by N, and £ is a scaling factor assuring that g(r) — > 1 
for r — > oo. Due to the finite integration range, all pairs of particles with a distance 



rij G [r — 1/2, r + 1/2] contribute to g(r). In Fig. 11 'b), the time dependence of g(r) is 
shown for a spherical interface and x — 0.305. One notices that all maxima decrease 
in time. This means that the local ordering decreases, particularly in the first 10 4 time 
steps. Afterwards, g(r) changes only slightly. Fig. [TTfc) compares g(r) for x = 0.153 
and x — 0.305 after 10 3 time steps and can be used to explain the differences of the 
snapshots in Fig. [To{a) and (b): the order is higher for a smaller concentration x- 
In the case of x = 0.305, the ellipsoids are attracted by capillary forces leading to 
a clustering of particles. This causes disorder, which manifests itself in smaller peak 
amplitudes of g(r). 

In conclusion, the particle coverage x aas a marginal influence on the duration 
of the flipping towards the interface. However, the particles rotate faster when the 
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Fig. 11: (a) Time evolution of the uniaxial order parameter 5 (defined in Eq. ([9j) 
the order decreases faster for the case of a sphere while the coverage fraction does not 
play an appreciable role, (b) The time dependence of the pair correlation function 
g{r) (defined in Eq. (10)) for particles adsorbed at a spherical interface (x = 0.305) is 
shown for different times, (c) Pair correlation function for the sph erical interface after 
10 3 time steps for different values of x corresponding to Fig. 10 'a) and (b). One can 
see that the order for the denser case is reduced due to capillary interactions which 
cause the particles to cluster causing a reduced long range order. 



interface is spherical. For increasing the spatial particle ordering is reduced since 
the average particle distance is sufficiently small to allow for attractive capillary forces 
which in turn cause particle clustering. 



3.3.2 Effects of particles on a droplet in shear flow 

In many industrial applications, the particle-covered droplets as discussed in the previ- 
ous section are not stationary or in equilibrium, but are instead subjected to external 
stresses or forces. The properties of the individual droplets are then of interest as 
their behaviour dictates that of, for example, an emulsion formed of these droplets. 
In this example, we consider monodisperse neutrally wetting spherical particles. The 
particles are adsorped to the droplet interface which is initially spherical. As before, 
we employ the particle coverage fraction the fraction of the initial interfacial area 
taken away by the presence of these nanoparticles. In the simulations, variation of 
X is achieved by varying the number of particles N while keeping their radius r p 
constant. External shear is realised by using Lees-Edwards boundary conditions [79], 
adapted for use with lattice Boltzmann [80 . As this shear is applied to the system, 
the droplet will start to deform. To quantify this deformation, for small to moderate 
shear rates 7, a dimensionless deformation parameter is used: D = (L — B) / \L + B) 
where L is the length and B is the breadth of the droplet. This parameter will be zero 
for a sphere and tend to unity for a very strongly elongated droplet. As the droplet 
loses its spherical symmetry, it will also start to exhibit an angle of its long axis with 
respect to the shear flow: the inclination angle 9d- 

Increasing the shear rate a droplet is subjected to will generally increase its defor- 
mation, as well as reduce its inclination angle from an initial angle of 45 degrees: the 
droplet is elongated and aligns with the shear flow. However, this is only valid as long 
as inertia can be neglected. When inertial forces are comparable to or stronger than 
viscous forces, the inclination angle can first increase beyond 45 degrees, before its 
eventual reduction. We now discuss the effect of increasing the particle coverage frac- 
tion x at constant shear rates 7. Representative snapshots of the droplets for various 
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Fig. 12: Effect of adding monodisperse and neutrally wetting spherical nanoparticles 
to the inferface of a droplet in shear flow. Representative snapshots of deformed 
particle-covered droplets are shown for 7 = 4.7 ■ 1CT 4 and (a) x = 0, (b) x = 0.27, 
and (c) x — 0-55 ([17], reproduced with permission of the Royal Society of Chemistry). 
The relative deformation of the droplet increases strongly for high coverage fraction, 
as shown in (d). In turn, the relative inclination angle decreases, indicating a better 
alignment with the shear flow, as shown in (e). 



X and fixed 7 = 0.47 • 10 -3 are presented in Fig. 12 a), (b), and (c). The particles are 
not homogeneously distributed over the droplet surface when shear is applied, due 
to an interplay between local curvature and shear velocities: high curvatures and low 
velocities are energetically favoured. 

Fig. [l2^d) shows the change in deformation of the droplet as the particle coverage 
is increased. The deformation has been rescaled as D* = (D — D°)/D°, where D° = 
D(x = 0). Due to a large decrease in free energy granted by the presence of particles 
at the interface, they are irreversibly adsorbed. When shear is applied, the particles 
are affected by this, and they start to move over the droplet surface, but cannot be 
swept away. As can be seen from the figure, high particle coverage fractions lead to a 
large increase in deformation of the droplet. The inertia of the particles plays a critical 
role here. The particles have to reverse direction to stay attached to the droplet, and 
massive particles strongly resist this change in velocity. This causes the interface to be 
dragged by the particles, increasing elongation from the tips. Fig. 12 [e) demonstrates 
the effect of particle coverage on the inclination angle of the droplet. Similar to the 
rescaling performed on the deformation, the inclination angle has been rescaled as 
®d = (@d - Qd)/Q% where 9° d = 9(x = 0). The inertial effects described above also 
cause the inclination angle to strongly decrease for large x- 
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4 Conclusions 

We have presented a non-exhaustive collection of simulation methods which are able 
to catch the relevant physical ingredients to study the dynamics of different kinds 
of fluid-fluid interfaces. The aim of this contribution is not to provide all strengths 
and weaknesses of these methods, but to demonstrate physical systems, where the 
selected methods are particularly suitable. Obviously, numerous alternative simula- 
tion approaches exist, but it is beyond the scope of this article to judge them by 
systematic quantitative benchmarks. 

The immersed boundary method offers a flexible and simple approach to model 
deformable particles immersed in suspending fluids, although it is known to be some- 
what problematic in the case of rigid boundaries, due to its explicit algorithm. As 
it is a front-tracking method, the interface configuration is directly known, which 
simplifies the force evaluation based on the interface deformation. Another advantage 
is that the constitutive properties of the interface can be controlled directly without 
tuning additional simulation parameters, and there is no need to remesh the fluid do- 
main. The IBM is especially suitable to simulate the dynamics of vesicles and blood 
components at small and intermediate volume fractions. 

The combination of the lattice Boltzmann method with the Shan-Chen multi- 
component model for the simulation of immiscible fluids and a lattice-based sur- 
factant model allows to coarse-grain the molecular interactions sufficiently well so 
as to reduce the computational effort substantially. Thus, large system sizes can be 
reached allowing to study for example macroscopic rheological properties and their 
dependence on microscopic interactions. Furthermore, the model does not require any 
information of the expected properties of a fluid mixture and the corresponding fluid- 
fluid interfaces because all interactions are included on a mesoscopic level providing a 
purely kinetic approach to the description of complex fluids. This allows, for example, 
to study the formation of cubic mcsophases such as the gyroid phase without adding 
any information on the gyroid itself a priori. However, a drawback of the model is 
that its phenomenological forces and parameters are often difficult to be related to 
experimentally available observables and time consuming parameter searches might 
be required to reveal interesting physical behaviour. 

When the multi-component lattice Boltzmann method is coupled to a molecular 
dynamics algorithm for the description of suspended (colloidal) particles, particle- 
laden interfaces can be studied at a level where not only hydrodynamics, but also 
individual particles and their interactions are resolved. At the same time, efficient 
scaling to allow larger domains to be simulated is still provided due to the locality of 
the algorithm. While this method has been proven to be particularly suitable to study 
particle stabilised emulsions, one of its drawbacks is the diffuse interface between 
fluids in the Shan-Chen LBM. The particles have to be substantially larger than these 
diffuse interfaces in order to be able to stabilise them. Thus, even comparably simple 
applications as the ones provided in this article require large lattices and therefore 
access to supercomputing resources. 
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